clc; clear; close all;
% 展示在不同耦合强度 k 下，x 与 y 的分岔行为，表现系统的周期、混沌等多种动力学状态

% 参数设置
mu = 3.5;
x0 = 0.1;
y0 = 0;
z0 = -1;

nk = 5000;                          % k 的采样点数
k_vals = linspace(0, 0.2, nk);      % k 扫描区间 [0, 0.2]

N_transient = 500;   % 瞬态迭代步数，消除初始影响
N_bif = 100;         % 分岔数据采集步数

% 初始化存储 x 与 y 分岔数据的单元数组
x_data = cell(1, nk);
y_data = cell(1, nk);

x_color = [186, 85, 211] / 255;
y_color = [255, 182, 193] / 255;

% 扫描每个 k 值并采集分岔数据
for i = 1:nk
    k = k_vals(i);  % 当前的耦合强度 k
    
    % 初始化状态 [x, y, z]
    state = [x0, y0, z0];
    
    % 瞬态迭代：进行 N_transient 步迭代，消除初始条件影响
    for n = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    
    % 分岔数据采集：再迭代 N_bif 步，记录每一步的 x 与 y 值
    x_temp = zeros(1, N_bif);
    y_temp = zeros(1, N_bif);
    for n = 1:N_bif
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        x_temp(n) = dx;
        y_temp(n) = dy;
    end
    x_data{i} = x_temp;
    y_data{i} = y_temp;
end

% 绘制分岔图
figure;
hold on;
% 遍历每个 k 值，将对应的 x 和 y 数据在 k 轴上以散点图形式绘出
for i = 1:nk
    k_val = k_vals(i);
    plot(k_val * ones(size(x_data{i})), x_data{i}, '.', 'Color', x_color, 'MarkerSize', 2);
    plot(k_val * ones(size(y_data{i})), y_data{i}, '.', 'Color', y_color, 'MarkerSize', 2);
end
xlabel('k');
ylabel('x, y');
title('x, y 随 k 的分岔图');
grid on;
